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Abstract. 

This paper describes a solution to the flexural analysis 
of a beam by digital oomputer. Subroutines were written 
which calculate shear, moment, slope and deflection at any 
specified point on a beam if the loading is known. Shear 
deflection may be included for many support conditions. 
Neither elastic supports nor column effects are considered, 
The subroutines solve problems of variable beam cross-section 
and loading by utilizing an extension of "McCaulay's Method", 
a generalized step function. 

The main program provides input-output facilities and 
also solves for indeterminate reactions on the beam, 

Problems and their solutions, showing the versatility 


of this program, are also inoluded, 
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1. Introduction. 


The routine solution to engineering problems by digital 
computer appears to be limited only by the scope of the 
available programs. A flexible program permitting the solu- 
tion of routine problems in flexural analysis of beams would 
seemingly fill a need in the field of digital computer ap- 
plications to Mechanics. The present work is devoted to such 
an idea. 

Five subroutines, "LOAD", "SHEAR", "MOMENT", "SLOPE", 
and "DEFLECT", were written which will calculate the load, 
shear, moment, slope, and deflection, respectively, at any 
indicated point on a beam if the indeterminate quantities 
are calculated elsewhere. The use of generalized step functions 
in this project allows for piecewise linear variation of 
distributed loading, bending compliance, =, and shear conm- 
pliance, 3. 

The main program, "BEAM3", first generates and then 
solves a set of simultaneous linear equations, calculating 
the indeterminate quantities of the beam, The above mention- 
ed subroutines are then utilized to calculate the remaining 
flexural quantities. This program also provides input-output 


facility. 
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2. Definition of the problem, 

The problem of structural beam analysis is one in 
which the external loading applied to the beam is well de- 
fined, as is the geometric configuration and the nature of 
the beam itself. Where information about such a beam is re- 
quired, this information will most likely fall into the 
following categories: 
. internal shear forces, V, and moments, M, 
. Slope,é. 
Deflection, y. 


. Redundant reactions, either forces, Ry, or momenta, 
RM 
Z° 


Wh FF 


"Flexural analysis" is used here as meaning the process by 
which this information 1s generated from the known data, 

In attacking this problem of flexural analysis, the 
intention was to limit the scope as little as possible. 
However, several customary simplifications and limitations 
were made, First, elastic supports and beam column effects 
were not considered. Limitations concerning the type and 
number of external loads and the geometry of the beam will 
be discussed later, Also tne basic relations which follow, 
inherently restrict the deflection of the beam to small values 
and the internal stresses to values below the elastic limit 
of the structural material. Shear deflection, which is 
generally neglected in a problem of this type, has been 
included in some problems. 

a. Bending deflection (y,). 

The Bernoulli-Euler-Navier equation governing the elastic 

curve of a beam due to bending, consistent with the sign 


convention shown in Appendix B igs: 


2 
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and we use the approximation, curvature ~ a“y/ax*. 


= curvature = (d‘y/dx”) [a*(aysaxy?] >” 


The following equations, derived from the above equation 
and shown in Timoshenko Ve Me with differences due to the chosen 
sign convention, describ? the flexure of a beam (See Appendix B 
for notation). 

w= a(x), 

v = {wax , 

M = /Vax , 

eC, - 21 - (4, ax 

yy) = #9; ax , 

b. Shear deflection (yo). 

The equation governing the shear deflection of a beam, 
also given in Timoshenko /5/, compatible with the stated 
sign convention becomes: 


dyo _ _ XV 


ax AG 
It is evident that: 


/ KV) 
ya = /\aG/ ax 


c. Total deflection (y). 
The equations used in this thesis are the result of the 


sum of the above shear and bending deflections. They are: 


W = q(x) , 


1 
Numbers so indicated refer to the bibliography, 


3 
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0 = es: dx - —~ + C_## 
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Y = S ax + Coxe 


Note that step-discontinuities may appear in all of the 
above equations, except in the equation for deflection. These 
discontinuities limit the methods that can be used in handling 
the eauations., 

Three methods that could be used to solve the equations 
are a pure analytical approach, a step-by-step numerical 
integration scheme, and the "curly bracket" method. Of 
these methods, the analytical approach is generally restrict- 
ed to a specific problem and is unsuited for adaptation to 
digital computer programming. The numerical integration scheme 
can be used in the computer solution of flexural analysis 
and may have been programmed; however, our search of literature 
has not uncovered any such program. The "curly bracket" 
method was chosen for use in this thesis because it was felt 
that it would result in a more exact solution requiring less 

*In these equations, ZF, 1 onde, , represent the sum of 

ve 


we point forces and moments, respect ly, applied to the 
eam, 


**C, and C. are constants of integration, 





7 


















ae oe 80h wh 





 —— (1 ome eee mee ee at ‘ 


ALD a ll lll ee a ly 
(gt Ci EF ee © Sd bene 
3 ee ge ee 
eee 1 aes oe or 
ee — te oe @ ee ee el 
fll Gm 
oC — ae -——_— ee a 
ea.) «& eh 
——_ OS i ee ee a ee ee Pome 
_ «= = = ««« —— =< =. omy 


‘OSte—<{—<—_{<—-_»— —“*§<—- —_— ea « — ——_ 












— Si. ie = re aN. ag 


ee)? io? oe 





computer time, particularly if results were required at only 


& few points. 





3. "Curly bracket" method 
The "curly bracket" method is based upon the properties 


of a convenient notation defined ag follows: 


fx-a}" = 0 if x<a 
= (x-a)" if x>a 


where n is an integer=0. 


Appendix A contains further properties of expressions contain- 
ing curly brackets. Tne notation is generally credited to 
McCaulay /l/, but similar, though less useful notations can 
be traced to A.Fdppl and St. Venant. Also, the ideas in- 
volved are related to the use of the Lavlace Transformation 
with the result that this method or approximationg to it are 
reinvented frequently /8/ /9/. 

Using this method allows one equation, valid over the 
entire length of the beam, to be expressed for each of the 
basic eauations of flexure. The particular value of this 
method is that it permits the integration of these equations 
to be carried out easily and has the virtue of adjusting the 
constants of integration most conveniently. 

In addition, for use with a digital comouter, the "switch- 
ing" property of the curly bracket can be easily handled in 
FORTRAN by use of an "IF" statement. 

A simple example should convince most readers of the 
convenience of this method. Examine the following problen, 
which consists of a uniform beam, simply supported at each 


end with a concentrated load at the center of the span. 
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stant variation of bending compliance, In this thesis 


1 

El’ 
the method has been further expanded to include beams with 
piecewise linear variation of both bending compliance and 


shear compliance s 
= AG: 





4, Implementation. 

Although the formulas shown in Appendix A could be ex- 
panded further, a piecewise linear variation of distributed 
loads, bending compliance, and shear compliance was consider- 
ed capable of providing sufficient accuracy and flexibility 
for the resulting program. 

Fig. 2(a) illustrates a distributed load on a beam 
whose origin is at x=0. Fig. 2(b) shows an element of the 
distributed load. 


y 
x 

Fig. 2(a) 

y Ff - 4 
Qa 
x 
a—s 

=n 

Fig. 2(b) 


The analytical expression for this element is: 


Q(x) = Q, Sx-A} © + “pe {xa} * _ Q, {x-B} © - so fx-B} + 


The same relationship is used for an element, f(x), of bend- 
ing compliance and an element, g(x), of shear compliance upon 
substituting f for Q and g for Q respectively. 

The basic equations of flexure can now be expressed in 


the same manner as they were utilized in programming. They 
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are, 


NOL 
W(x) = > Q, (x) (1) 
i=l 


(x) { ( = ° a ae ° 
V(x = W(x)dx + F Xx-X os R X-X 
P i=1 y,t is {=i yt r,3 


(2) 


where Fy are upward, known forces at x=Xp and RY are up- 
1 


ot i 
ward, unknown reactions at X=X,, 
i 
Xx i 5 — e 
M(x) = § Vix)dx + M x-X + RM Xx—X 
O i=1 7 m,$ 1=1 241 roy} 


(3) 
where Mo are clockwise, known moments at X=X,, and RM, are 
i i 


clockwise, unknown reactive moments at XX 
i 


g dy x NOI NOK 
= oe. = M(x f,(x)adx -|V(x x)) + C 
ax fe | pa 1 am =; a . 
(4) 
where Cy ig the constant of integration. 
Appendix A indicates how the product and the integration of 
the product of the curly bracket symbols are formed, 
XK -Xe- NOL X- NOK 
y = M(5) >. f,(§ )} as ax - Ux) > g(x)ax + OC.x +O 
= re 1 2 
O i=l] IG, {=1 
(5) 
where ©, is the constant of integration. 
Equations (1) through (5) when properly coded in FORTRAN 


“The numbers, NOL, NOF, NOR, NORM, NOI and NOK, appear- 


ing above the summation signs are program input parameters 
indicating the number of terms for each problem. 


10 





form the core of the five basic subroutines: LOAD, SHEAR 

MOMENT, SLOPE and DEFLECT. It must be noted before discuss-— 

ing the subroutines that equations (2) through (5) involve 

the unknown values C,, C,, R, , and RM_ , where there may be 
1 2 y; Zs 

as many as ten values each of Ry and RM,. The calculation of 

these values is done in the main program BEAM3 and will be 


discussed in section 6, 


Bl 





5. The basic subroutines, 


Subroutines LOAD, SHEAR and MOMENT, 

The above named subroutines, as their names imply, are 
employed to calculate the total distributed load, shear and 
moment, respectively, at any desired point on the beam. The 
expanded forms of equations (1), (2) and (3) were used in 
writing these subroutines. Utilizing the distance, X4y, 48 
the entering argument, the subroutines select and sum the 
appropriate terms in the corresponding equation by inspect- 
ing all loading terms and discarding those in which the value 
within the curly bracket is negative, Actually the load, 
shear and moment are calculated a small increment, E = 1078, 
to the right of X; in order to obtain the right side value 
of any step-discontinuities occuring at the point, Xs (See 


flow charts on pages 13, 14 and 15). 


Subroutine SLOPE. 

This subroutine is equation (4), excluding the constant 
of integration, C,, in FORTRAN language (See flow chart on 
page 17). Note that the equation can be expressed as a 
summation of terms, which are constants times the multiplica- 
tion of two step-functions of the form C fx-a} " fx-b?”, or the 
summation of the integral of the same type of terms, where n 
equals either 0 or l. 

Function subroutines, FTON and FT1LN (see Appendix A), 
will calculate the value of fx-a}° {x-b}™ and fx—a} + {x-b} 2, 
respectively. Function subroutines, ENTON and ENTIN (See 
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Flow Chart of LOAD 









SUBROUTINE | 
LOAD(X,W) | 


pane 


ePsx=107 © 





WW # ly 


he a 








Note: Numbers preceded by a number sign, #, refer to 
the statement number in the listing of LOAD, 
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Flow Chart of SHEAR 


SUBROUTINE | 


SHEAR(X, V) 


EPSX#107° 
- + on XeEPSX + 
~XFO(N) 


VV=FY(N) 


NeNel 


a ae . : 
ee x a(t) (2) 
“i | 
ew Sa 


N-NOF 


0] 
~<X+EPSX-B(I) > = = 


we 6 | 





VaVeVV 





I-11 | 





RETURN 
1 - 


Note: Numbers voreceded by a number align, #, refer to 
the etatement number in the listing of SHEAR, 
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Flow Chart of MOMENT 






SUBROUTINE 
MOMENT (X,AM) 


EPSx=107° 






+ EPSX=XK (K ) > 


AMeAM+AMM 


{= Xe EPSX = 
~XFO(N) 
AMM= #10 







AMeAMSAMK | 











| AM#AM+AMM | 


ls 


[ mea | 





KA-NORM a 
, | 
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RETURN 


Note: Numbers precedec by @ number aign, #, refer *- th 
statement number in the listing of M°CM=NT, 
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Appendix A), will calculate the first integral of these 

functions. Equation (4) can now be expressed as the summation 

of constants multiplied by the appropriate function subroutines, 
SLOPE calculates all combinations of M(x) and f(x) by 

stepping through all the f(x) terms and varying all the terms 

in the moment equation at each step. It then steps through 

all the g(x) terms varying all the terms in the equation for 

shear at each step. The sum of these terms is the value of 

the slope, minus the integration constant, 0), at the entered 


distance, x. 


Subroutine DEFLECT, 

At this voint the author wishes to point out the similar-— 
ity of the equations for slope and deflection. Since the 
integration constants, Cy and Co, are again excluded from the 
subroutine, and the function subroutines, DITON and DITIN 
(See Appendix A), calculate the second integral of the mul- 
tiple step-functions, subroutine DEFLECT is identical to 
SLOPE except for substituting ENTON for FTON, ENTIN for FTIN, 
DITON for ENTON and DITIN for ENTIN. The result is a sub- 
routine that will calculate the deflection, minus the con- 


stants of integration, at any required point on the bean. 
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6. Program BEAM3. 

The main program, BEAM3, can be divided into four dis- 
tinct steps. These steps are: reading the input data, cal- 
culation of indeterminate quantities, calculation of the 
flexural quantities, and the output of the solution, 

The program first reads the input data which includes 
the length of the beam, all external loading, the location of 
all reactions, the flexural properties of the beam, and a group 
of parameters giving the number of each type of load or re- 
action, It also provides for reading a predetermined value 
of deflection, PDEFT, at each reaction and a predetermined 
value of slope, PSLP, at each moment reaction, 

The program then proceeds to solve for the indeterminate 
quantities of the beam by utilizing the two equations of 
static equilibrium plus equations obtained from the known 
boundary condition at each reaction or reactive moment. The 
first two equations are simply the summation of the forces 
on the beam equals zero and the summation of the moments 
about the right end of the beam equals zero, The remainder 
of the equations are obtained from the fact that the slope 
at each reactive moment or the deflection at each reaction 
is known. Using this fact in the appropriate one of either 
equation (4) or (5) results in a set of linear equations for 
the indeterminate quantities. 

The equations thus obtained are generated in a matrix 
form that is compatible with subroutine GAUSS2 (See Appendix £). 


This subroutine is then called to solve the equations simulta- 
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neously, If the matrix is singular, GAUSS2 exits to the call- 
ing program which in turn prints "MATRIX SINGULAR" and then 
stops. If not, the subroutine solves the equations and re- 
turns the calculated values to the main program, The most 
serious limitations of this program are imposed by the gen- 
eration or solution of the simultaneous equations, therefore 
it is appropriate to include and explain several of them at 
this point. 

First, at least one reaction, even if it is later forced 
to equal zero, must be included with the input, or a singular 
matrix will be generated, 

Including shear deflection in the solution must also be 
greatly limited at this point. Since in general the shear 
includes step discontinuities, they will also appear in the 
slope if shear deflection is included, The program generates 
an equation for slove at each moment reaction. If by chance, 
either a reaction or point force occurs at the same point as 
a moment reaction, the slope at this point is double valued, 
The program is unable to distinguish which of the two values 
is required except when it occurs at the extreme ends of the 
beam, For example a beam "built in" at both ends is accept- 
able. If the beam is "built in" at only one end, it must be 
solved in a manner such that the “built in" is at the origin 
of the bean, 

The program then divides the beam into a number of 


increments, NOP-1,* and calls the subroutines, LOAD, SHEAR, 


“NOP is another parameter read in at the start of BEAM3, 
We, 
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MOMENT, SLOPE and DEFLECT to calculate the flexural quantities 
at each of these increments, 

The output consists of printing the distance from the 
left end of the beam to the point in question and the values 
of the load, shear, moment, slope and deflection at that 
point. The values of all of the indeterminate quantities are 


also printed out. 
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Fig. 7 Simplified Flow Chart of BEAM3 
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7. Testing of Program, 

The testing of the program and the subroutines could, 
in theory, continue indefinitely because of the unlimited 
number of beam configurations, The author has attempted to 
show (See Appendix D) solutions to test problems of simple 
configurations and constant bending and shear compliance to 
compare with classical solutions. fThen a statically deter- 
minate problem with a varying bending compliance was solved 
and compared to a solution obtained from a numerical in- 
tegration scheme, The most rigid test for the program was 
in constructing a problem for a beam of unlikely configuration 
(See Problems 6 and 7 of Appendix D), and solving the flexural 
analysis of this beam from both ends. Note that the only 
significant difference between the two solutions is the opposite 


signs for shear and slope, which would of course, be expected, 
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8, Conclusions and Recommendations, 

The thesis shows a method of including piecewise linear 
variation of both bending compliance and shear compliance in 
the flexural analysis of beams by extending "McCaulay's Method", 
Adapting this idea to the digital computer resulted ina 
program capable of analyzing a very general type of beam, 
However, had time permitted, several changes would have been 
made to the program, First, internal monitoring of the input 
data would be an asset to any user. Also, using the same 
general method of attack, this program would have been ex- 
panded to include bending compliance and shear compliance of 
a piecewise parabolic nature. 

Other additions to the program, much larger in scope, 
would be including elastic supports and/or beam column effect. 

It would be interesting to compare the results of this 
program to those obtained from a program utilizing the numer- 
ical integration method, Also, a more exhaustive study of the 
program could be made to determine accuracy of, and time 


requirement for, solutions. 
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APPENDIX A 
"Curly bracket" method, definitions and formulas 

We define: 

ky feeacy f = Gea) eh x>a 

= 0 ; x<a 
where n is an integer=20 

From this definition, it is possible to establish the 
following formulas, Although it would be possible to express 
these results in terms of curly brackets, in most cases the 
resulting expressions would be rather complicated. The forms 
exhibited below have the advantage of being most readily 


conformable to the decision commands available in FORTRAN, 





+1 
~ n (x-a)® 
é. j jx-a) ax " nel , we>e 
re) 
= 0 , Otherwise 
36 fx-aj® {x-b} © ~ (x-b)" ; xXx>a and x>b 
= 0 : Otherwise 
r 1 
4. §x-a} *{x-b}” : (x25) = (b-a)(x-b)" , x>a 
and x>b 
~ 0 : Otherwise 
n+l n+1 
-X , (x-b) (a-b) 
O n a (x-b) : 
5. \ §x-a} {x-b} ax ary ee , xX>a>db 
nel 
2 X-2 , x>b>a 
nel ee 
= 0 , Otherwise 
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Gas n (pa) Geme 


* aL n 
6. \ xa} * fx--b3 Pax 


n+e2 nel 
‘ toi) lle Meee: 
(n+1) (n+2) oss 
nt¢2 nel 
; - (x-b) > (bap ean ee 
ne2 n¢l 
= 0 , Otherwise 
eee - = 
te \ \, 96 -af $6 -b} ag ax 
2 
(xb) (Fenian ae 
n¢1) (n+2) (n+l) 
ne2 
S MeeD) a. X>a>b 
(n+1) (n+2) 
nt2 
_ (xb) 
(n+1) (n+2) x >b>a 
= 0 : Otherwise 


8. { ‘i ‘é -s}* fe -bp ae, dx 


nt3 








(x-b) (b-a) (x-b)"** 
(n+2) (n+3) * (n+1) (n+2) 
PEG 2(a-b) "* 
= (227) nee) i (n+1) (n+2)(n+#3) , x>a>db 
3 ne2 
5) mae (b-a) (x-b) 
~  (n+2) (n+3) y (n+l) (n+2) : eg 
= 0 , Otherwise 


Function Subroutines FTON, FTIN, ENTON, ENTIN, DITON, DITIN 
The above named function subroutines were written and 
utilized to evaluate each of the preceeding formulas numbered 


3 through 8. The values x, a, b andn were used as entering 


26 





arguments to these functions, The function then became the 


variable as follows: 


FTON =  $x-at® §x-b}™ 

FTIN. = {xa} {x-bz” 

ENTON = i, fx—e3 ° §x-b} Max 

ENTIN = A fx-a} * Gx-d} Max 

ETON = 3 (6-03° §e-ppPag ax 
DITIN = > " e-a}> {6 -vpPag ax 


As examples, flow charts for FTON and ENTON are shown in 
Figs. 8 and 9. 
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Fig. 8 Flow Chart of FTON 


FUNCTION 
FTON(X,A,B,N) 


FTON=0 








Fig. 9 Flow Chart of ENTON 


FUNCTION 
ENTON (X,A,B,N) 





4 X-A ; 





aon FNeL oO OFN#1 





RETURN 
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APPENDIX B 
Sign Convention and Notation 
Sign Convention 


The following sign convention has been employed: 


’.. y = deflection 


slope = ay 


ax 


M (| a ) M = bending moment 
V {| i V = shear force 
i St 


Note: 1. Point forces, F,, and point reactions, R 
positive in the positive y direction, 


\ 
D 
al 
D 
i] 


distributed loads 


y» are 


2, Point moments, M,, and point reactive moments, 
RMo, are positive in a clockwise sense, 


Notation: 
THESIS PROGRAM 
C Length of beam 
W W Total distributed load at a point 
_V V Shear 
M AM Moment , 
S, Yl Slope 
y Y Deflection 
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f(x) 


@(x) 


RY 
AMZ 
RAMZ 


AQ 


GA 
GB 
XFO 
XR 
XM 
RXM 
AG 
BG 
AE 
BE 
A 

B 
NOL 


Point load or force 

Point reaction 

Point moment 

Point reactive moment 
Element of distributed load 
Distributed load (left end) 
Distributed load (right end) 


Element of bending pte, in 
= ) 
Bending compliance EI 
(left end of element) 
Bending compliance 
(right en¢e of eae ha 
Element of shear comp asgce 


(E) 
Shear compliance 

(left end of element) 

Shear compliance 


(right end of element) 
Distance to Fy 


Distance 
Distance 
Distance 
Distance 
Distance 
Distance 
Distance 
Distance 


Distance 


to 
to 
to 
to 
to 
to 
to 
to 


to 


R 


td 


M 


aan 


ro) 
p 


p oF 


ED 
Qo 
Qp 


Number of distributed loads (25)? 


Numbers following explanation of terms NOL--NOK indicate 
the arbitarily chosen maximum value of each for this progran, 
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NOF 
NOR 
NOM 
NORM 
NOI 
NOK 
PDEFT 
PSLP 
NOP 


NOF 
NOR 
NOM 
NORM 
NOI 

N OK 
PDEFT 
PSLP 
NOP 
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Number of point forces (100) 
Number of point reactions (10) 
Number of point mom3nts (50) 


Number of point reactive 
moments (10) 

Number of bending compliance 
elements (25) 

Number of shear compliance 
elements (25) 

Predetermined deflection 

at a reaction 

Predetermined slope at a 
reactive moment 

Number of points at which 
the flexural quantities will 
be calculated, 





APPENDIX CG 
General Information Concerning Use of Program 

The composite program, BEAM3 plus subroutines, was 
written in FORTRAN-60, compiled and tested on a Oontrol 
Data Corporation 1604 digital computer. It consists of 412 
FORTRAN statements requiring about 580 cards, The computer 
storage requirement is about 12,500 cells, however, this 
number can be reduced by changing the dimension statements 
which would of course, reduce either the maximum number of 
external loads, elements of shear compliance or bending compli- 
ance, or number of points at which flexural quantities could 
be calculated, 

The FORTRAN statements used in this program are: GO TO, 
computed GO TO, IF, STOP, DO, CONTINUE, FORMAT, READ, PRINT, 
WRITE OUTPUT TAPE, FUNCTION, SUBROUTINE, RETURN, CALL, DIMENSION, 
COMMON, and END. 

There are no conversion constants incorporated in this 
program. As such all input must be in a consistant set of 


units, 


"FORMAT" for Input Data 

The input data for a problem is read by the program 
with a FORMAT of I4 for all "fixed point" variables and F20,0 
for all floating >voint variables, The arrangement of the 


data on the input cards is shown as follcws; 


BI. 
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Group 
no. 


1 


—- WwW WW 


oOo oO NN BD WA 


No. 


of cards 


in group 


L 
L 
NOL 
NOF 
NOR 
NOM 


NORM 


NOT 
NOK 


Fields 
F20.0 
814 
4F20 .0 
2F20.0 
3F20.0 
2F20.0 
3F20.0 
4F20.0 
4F20 0 


Entries 
C 


NOP, NOL, NOF, NOR 
NOM, NORM, NOI, NOK 
AB, Ag, wi: 

Xf, Fy 

Xp, 0.0, PDEFT 


Xn» Mz 


KX, 0505, Pole 


m? 
he, Ber Eg, EB 
Ap, Bos Ga» G, 


As an example the input data cards for test problem #1 


in Appendix D was ag follows: 


CARD 
CARD 
CARD 
CARD 
CARD 
OARD 
CARD 
CARD 


100 


ou lye 0” B2arO 2 


0. 
0. 
100 
0. 
100 
QO. 


100 
0. 
0. 
0. 
0. 


100, .000000001 
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0 


-~ 500. 


oO 
oO 
0. 
0. 


~ 500. 


»000000001 








aT 


tt e+ Tete 





‘tet 


wee . 
aa 


APPENDIX D 


Test problems and solutions 


Test problem l 
This problem consists of a beam of uniform cross-section, 


"uilt-in" at both ends, and loaded with a uniformly dis- 


tributed weight of 500#/in. (See Fig. 10) 


——— les iE = 
BI 2 #in.* 
K 

AG * 





Test problem 2 
This problem is identical to problem 1 except that the 


slope at either end of the beam is forced to values other 


than zero, 


ej =o 
left end age 


fe} = F 2 
© sight end a 
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SOLUTION TO TEST PROBLEM 1 


SHEAR MOMENT SLOPE DEFLECTION 


LOAD 


DISTANCE 
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SOLUTION TO TEST PROBLEM 


SLCPE DEFLECTION 


LOAD SHEAR MOMENT 


“DISTANCE 
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Test problem 3 
This problem consists of a beam of uniform cross-section, 


"built-in" at one end and simply supported at the other, and 


loaded with a uniformly distributed weight. (See Fig. 11) 


| Qe ; 1 ; 
. (> O#/in ez = -1332x10 72, 
a ee : 
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Test problem 4 
This problem is identical to problem number 3 except 


shear deflection is included in the solution. 


K ~6 
iG = ot 59 x 10 1/#in. 
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SOLUTION TO TEST PROBEEM™ 3 


SLOPE DEFLECTION 


LOAD SHEAR MOMENT 


DISTANCE 
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SOLUTION TO TEST PROBLEM 4 


SHEAR MOMENT SLOPE DEFLECTION 


LOAD 


DISTANCE 
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Test problem 5 

This problem consists of a homogeneous circular shaft, 
three inches in diameter at the center and tapered-to one 
and a half inches in diameter at both ends, The total weight 
of the shaft is 10047 and a uniformly distributed load total- 
ing 200# is placed on the central section, The shaft is 
supported by a uniformly distributed reaction at each end. 
Points A and B are points of zero deflection, This problem 


was approximated for the program as shown in Fig. 12. 


8.33 Fin. 


[-—>| Yin, —— 


+A +¢ B+ 


=. i61n. 2 Fan. 1b in. da 


Fig. le 

Weight of shaft: 

o" — 8" -0,.5315#/in. 

8B" —~ 24" -0.5315 to -2.08l1#/in, (linear) 

24" — 36" ~2,.08l#/in. 

Bending compliance (=| ; 

or — 9" 0.1420x10_, 1/#-in* 6 5 

ge" . 12" 0.1420x10" - to 0.0583x107 1/#-in fm 
Tas— 16" 0.0583x107 7 to 0,0281x107 1/#-in linear) 
16" — 20" 0.0281xl07? to 0.0152x107 1/#-ing (linear) 
20%" 24% 0,.0152x107~° to 0,0089x107 1/#-iné (linear) 
2u" _ 36" 0.0089x10° 1/#-in2 


K 
Shear compliance (=a) mPO 
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SOLUTION TO TEST PROBLEM 5 


SHEAR MOMENT SLOPE DEFLECTION 
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DISTANCE 
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Test problem 6 
This hyperstatic problem consists of a beam with a very 
unlikely configuration as shown in Fig. 13. The bending com- 


pliance varies as shown in Fig, 14, Shear deflection has been 


K 
omitted, 1.¢., AG = Q. 
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DISTANCE Alone BEAM Fir, 


Fig, 14 


Test problem 7 
This problem is identical to problem number 6 except the 


beam was turned end for end, 
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SOLUTION: TO TEST PROBLEM 6 


SHEAR MOMENT SLOPE DEFLECTION 


LOAD 


DISTANCE 
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SOLUTION TO TEST PROBLEM 7 


LOAD SHEAR MOMENT SLOPE DEFLECTION 


DISTANCE 
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APPENDIX E 
Subroutine GAUSS2 

This routine is a modified form of the GAUSS2 used by 
C. B. Bailey /7/ in his program, LINEQN, for solving sinmul- 
taneous linear equations by Gaussian elimination, Bailey 
was solving the matrix equation, Ax#B, with a possibility 
of 50 vectors B for each matrix A, but this program requires 
only one vector B with each matrix A. Also since the maximum 
number of equations in BEAM3 is 22, the size of matrix A was 
reduced from 100x101 to 22x23. 

At each step of the elimination, the value of the diagonal 
element is compared to a defined "zero", If it is smaller 
than this quantity, the matrix is considered singular and 
an error output is returned to the calling program, This 


is the "matrix singular" test used by BEAM3, 
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APPENDIX F 


Listing of program and subroutines 


Name Page 
BEAM3 48 
LOAD 54 
SHEAR Ne, 
MOMENT 56 
SLOPE 58 
DEFLECT 61 
FTON and FTIN 64 
ENTON and ENTIN 65 
DITON and DITIN 66 
GAUSS2 67 
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PROGRAM BEAM3 

DIMENSION WC 1001) 4V(1001),4M(1001),Y1(1001),¥(1001), XX( 1001), 
1 A(25)_¢B(25) sQA(25),Q8(25) sXFO(100),FY(100)sXRE10)sRY(10)s 
2XM(50) »AMZ(50) RXM(10) »RAMZ( 10), AE( 25) ,BE (25), £A(25) »EB(25)4_ 
3AA(22423)yXA(22) » PDEFT (1 0) PSLP(10) 
kh sITITLE( 10) +GA(25)yGB(25) AG(25) ,86 (25) 

COMMON AyB,QAyQBy XFO, FY, XRy RV yXMy AMZ, RXMzRAMZy AEs BEy EAy EB yNOLY 
1 NOR,NOF,NOM,NORM, NOI ,NOK,GA;,GB, AG,BG 

READ MeeTON.C 


70 


71 


50 
2 
51 
52 
73 
Sys. 
54 
Th 
515 
56 
57 
58 
59 
60 
61 
62 


‘NOQ] 


FORMAT (F20.0) 

READ 71,NOP,NOL,NOF, NOR, NOM,NORM, NOI y NOK 

FORMAT(814) 

IF(NOL) 51,51,50 

READ 72,(A(1),B8(1),QA(I) »QB(1),1=1,NOL) 

FORMAT(4F20.0) 

IF(NOF) 53,53,52 

READ 73,( XFO(N) »FY(N) »N=1,NOF) 

FORMAT (2F20.0) 

IF (NOR) 55,55,54 

READ 74,(XR(M),RY(M) »POEFT(M) ,M=1,NOR) 

FORMAT (3F20.0) 

IF (NOM)57,57,56 

READ 73,(XM(K),AMZ(K) ,K=1 NOM) 

IF (NORM) 59,59,58 

REAN 74, (RXM(KA) »RAMZ( KA)» PSLP(KA) »KA=1, NORM) 
IF (NOI)61261,60 

READ 72, (AEC IA), BE(IA) EAC IA) »EB( IA), 1A=3,NOI) 
IF (NOK) 1,1, 62 
READ 72, (AG(IB) »BG( 1B), GA( IB) »GB( IB), I18=1,NOK) 

NORT = NOR+1 

a0. 

G2— 70. 

NORM2 =NORM + 2 

NORM3 =NORM + 3 

NOQ = NOR +NORM 

NOQ + 1 

NOQ + 2 


NOQ2 


48 





55 


WS 
ke 


U8 


Nows ="NOQ + 3 
EPSX=10.##(-8) 
DO 2 K=1,NCR 


ARA(1,K) = 1. 

DO 4 K= NORI,NOQ2 
AA(1,K) = C. 
AAA=C. 


IF (NOF)45,45,35 
DO & N=1,NCF 
AAA= AAA + FY(N) 
ITFCNOL)4E ,48,4S 
CO 5 f=1,NCL 
AAA = AAA + ((QA(T)+QB(1))/2.)% 
AA( 1,NCQ3)=-AAA 
O00 6 K =1,NCR 
AA(2,K) = C-XR(K) 
IF(NCRM)9,9,7 
DO @ K=NOR1,NO0Q 


8 AA(2,K)=1. 


11 


DO 1C K=NCQ1,NGQ2 

AA(2,K) = 0. 

CALL MOMENT (CysAAA) 

AA(2,NOQ5) =—-AAA 

IF (NORM)17,17,11 

DO 16 J=3,NORM2 

DO 12 K=1,NCR 

AA(J,K)=0. 

D0 12 TA=1,NOI 
AAA=EA( TA) #®ENTON(RXM(J-2), AEC TA) sXR(K) 91) 


(BCI )-ACT)) 


- EBCTIA) #ENTON(RXM(J-2), 


‘TBECIA) »XR(K),17) +CCEBCIA)“EACIA)) /(BECIA)-AE( IA) )) © (EATIN(RXM(J-2) 
2,AE(IA) oXR(K) 91 )-ENTIN(RXM(J-2),BE( IA) »XR(K), 1)) 


12 


900 


902 
$C} 


AA(SeK) = AAC J,K) + AAA 
IF (NCK) 13,13,900 

DO 91C IB=1,NOK 

IF (NORM-—1)902,901 ,SC2 

IF ( J-NORM2)901,904,901 


X=RXM(J-2) -EPSX 
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AAA = GA(IB)# FTON( X sAGT1S),XR(K),0) = COILS = F Von. 


1 »BG(18),XR(K) 20) + ((GBCIB)-GA(IB))/(EG(IBI-AG(IB)))« ( FTIN( 
2 x »AG(IB) sXR(K),O0)— FTIN( xX BG CIG 1s MRK eeO Is 
GO TO 910 
904 X=RXM(J-2) +EPSX : 
AAA = GA(IB)# FTON( xX »AG(IB) +XR(K),0) -— GB(IB)* FTON(X 
1 »EGCIB),XR(K) 0) + ((GBCIB)-GA(IB))/(BG(IB)-AG(IB)))* ( FTIN( 
a xX MNGUMENRREK);»Ol= FTING x BG (1B), RK ne om 


910 AA(J,K) =AA(JsK) -AAA 
12. CONTINUE 
IF ( NOC-NOR1) 155,135,135 
135 CO 15 K=NOR1,NCQ 
AA(J,K)= 0. 
DO 14 IA=1,NOI 
AAA= EA( IA) #ENTON(RXM(J-2) pAE( IA) »RXM(K-NOR),C) -EBCIA)® ENTCN 
T (RXM(J-2),BECIA) »RXM(K-NOR),0) + ((CEBCIA)-EA(IA)) /( BEC IA)-AE(IA)) 
2) #(ENTIN(RXM(J-2),AE( IA), RXM(K=NOR),0) -— ENTIN(RXM(J-2),BE( IA), 
3 RXM(K=NOR) ,O)) 
14 AA(J,K) = AAC JK) + AAA 
15 CONTINUE 
155 AA(J,NCQ1)=1. 
AA(JS,NOC2)= 0. 
~ ITF (J-NORM2) 156,157,156 
156 X=RXM(J-2) - EPSX 
CALL SLOPE(X,AAA) 
GO TO 16 
157 X=RXM(J-2) + EPSX 
CALL SLOPE(X,AAA) 
16 AA(J,NOG3) =-AAA + PSLP(J-2) 
17 DO 2C J=NORM3,NOG2 
00 1€ K=1,NOR 
AA(J,K) =0. 
DO 181 IA=1,NOI 
AAA = EA(IA)® DI TON( XR(J-NORM2),AE(IA) sXR(K)_1) -ERC(IA)*CITCN 
1 (XR(J-NORM2) »BEC IA)» XR(K)_1) +((EBCIA)-EACIA))/ (BECIA)-AECIA))) 
2 #(OLTIN(XR(J-NORM2),AE( IA) »XR(K)91) -— DITIN(XR(J-NORM2),BE( IA); 
3 XR(K)51)) 
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181 AA(JsK) = AA(J,K) + AAA 
IF (NOK)18, 18,920 
920 DO 930 IB=1,NOK 
AAA =GA(1B)#® ENTON(XR(J-NORM2),AG(IB) »XR(K);0) - GBCIB)#ENTON 
1 (XR (J-NORM2),BG (1B) sXR(K),0) #((G8(1B)-GA(TB))/ (BGC 18)-AGC1B))) 
2 *(ENTIN(XR(J-NORM2),AG(IB) »XR(K),0) -ENTIN(XR(J-NORM2)yBG(I8) » 
3 XR(K),0)) 
930 AA(J,K) =AA(J,K) -AAA 
18 CONTINUE 
IF(NOQ-NOR 1) 195,185,185 
185 DO 19 K=NOR1,NOQ 
AA(JyK) = 0. 
00:19 IA=1,NOI 
AAA = EA(IA)* DITON(XR( J-NORM2) »AEC(IA)»RXM(K=NOR),0) - EB(IA)# 
DITON( XR( J-NORM2) »BE( IA) »RXM(K-NOR) ,O) #((EBCIA)-EA(IA)) / (BEC(IA 
2 y= AE(IA))) #® (DITIN(XR (J-NORM2),AE(IA) »RXM(K-NOR),0) - DITIN 
3 (XR(J-NORM2),BE(1A) +RXM(K=NOR) 70)) 
19 AA(JyK) = AAC J,K) + AAA 
195 AA(J,NOQ1) =XR(J-NORM2) 
AAC J,NOQ2) = 1. 
CALL DEFLTCT (XR(J-NORM2), AAA) 
20 AA(J,;NOQ3) =-AAA+ PDEFT (J-NORM2) 
~ CALL GAUSS2 (NOQ2,.00000001;AAy XA,K1) 
GO TO (22:21),K1 
21 WRITE OUTPUT TAPE 4,77 
77 FORMAT (16H MATRIX SINGULAR) 


STOP 
22 DO 23 K=1,NOR 
23° RY(K) = XACK) 

IF (NORM) 26,26,24 


24 DO 25 K=NOR1,NOQ 
25 RAMZ(K=NOR) = XAC(K) 
26 Cl = XA(NOQ1) 
C2 = XA(NOQ2) 
X=0. 
NOPE =NOP-1 
FNOP = NOP-1 


pul 





DO 27 J=1,NCPE 
CALL LOAD (X,W(J)) 
CALL SHEAR (X,VCJ)) 
CALL MOMENT (X,AM(J)) 
Xl =X+EPSX 
SAUL SPOPECK1,¥ 10 J)) 
YI(JI=¥1( 5) +Cl 
CALL DEFLECT (Xs ¥f{J)) 
Y(J)=V (J) 40 1#X+C2 
xq) = ox 
27 X = C/FNOP +x 
XX (NOP) =x 
X=X-3.#EPSX 
CALL LOAD(X,WCNOP )) 
CALL SHEAR(X,V(NOP)) 
CALL MOMENT(X,AMCNOP) ) 
CALL SLOPE(X,YI(NCP)) 
YI(NOP)=Y1(NOP)4C 1 
CALL CEFLECT(X,Y(NCP)) 
Y(NOP)=Y(NOP)+C1#X 4C2 
WRITE OUTPUT TAPE 4,76 
76 FORMAT ( 3Xy1HJs7X+4HLOADs 10Xy SHSHEAR » 10X, 6HMCMENT 8X, 5HSLOPE, 
~ 1 9X, 1CHDEFLECTION, 8X, 1HX) | 
WRITE OUTPUT TAPE Uy7Ee(JaW(S) VES) 9AM(S) YC 5) VCS) 9 XX0S) del 
NOP) 
78 FORMAT(14,6E15.8) 
IF (NCR)282,282,260 
280 PRINT 281,(M,RY(M) »M=1,NOR) 
281 FORMAT ( 3HRY(,12,2H)=,E20.8) 
282 IF (NORM) 288,288,282 
283 PRINT 284,(KA,RAMZ(KA),KA=1,NCRM) 
284 FORMAT (SHRAMZ(,12,2H)=,£20.8) 
288 PRINT 286,C1 
286 FORMAT(5X,3HC1=,£20.8) 
PRINT 287,C2 
287 FORMAT (5X, 3HC2=,E 20.8) 
285 CONTINUE 
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SUBROLTINE LOAD(X »W) 

DIMENSION A(25) +B (25),QA( 25) »QB(25) ,XFOl 100) ,FY(100),xXR(10) »RY(1C) 
19XM(5C),AMZ(50) 2XM( 10), RAMZ( 10) ,AE(25) »BE(25) »EAL 25) yEB(25) 
2+GA(2£)+GB(25),AG(25),B8G6(25) 

COMMON Ay ByQAyQB, XFOs, FY, XR RY XMpAMZ yo RXMy RAMZ AE, BEy EAy EByNOLy 
1 NOR, NOF,NOM,NORM,NCI,NOK,GA,GB,AG, BG 

EPSX=10.%*(-8) 

W=0. 

D0 5 I#1,NOL 

IF (X+EPSX-B(1)) 3, 3,5 

IF (X+EPSX-A(1) 15, 4 yd 

WW=QA(T) #((QB( I )—GQACT) )# (XACT) (BC T)-ACT)) 

WeW+Wh 

CONTINUE 

RETURN 

END 
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SUBROLTINE SHEAR(X,V) 

DIMENSION A(25),8(25) ,QA(25),QB(25) »XFC( 100) ,FY(100),XR(1C),RY(10) 
1,XM(5C) ,AMZ(50) sRXM( 10), RAMZ(10) ,AE(25) ,BE (25), EAC 25), EB(25) 
2,GA(25),6B(25),AG (25) ,8G6(25) 

COMMON A, B,QA,QBy XFO, FY s XRyRYyXMyAMZyRXMyRAMZyAE VBE, EAVEByNOL 
1 NOR,NOF,NOM,NORM,NCI,NOK,GA,GB,AG, 86 

EPSX=10.2#(-8) 

V=0. 

IF(NOL) 7,791 

ieee 6 [=1,,NOL 

IF (X-A(I) 16,252 

2 IF(X+EPSX-B(1))4,4,2 
3 VV= QA(I) #(X=$A(T) 240 (OBC T)-QACT)) #C(X-ACL) )##2))/02.8( B01) -AUT))) 
T-QB(1)#(X-B(1T)) -( (QB (1 )-QACT)) #((X-BC 1) ) ##2))/(02. (BC TI-A(L))) 
GO TO 5 

VV=4QA(T) #(X-ACT) 24° (QB( IT )-QACT) 9 # ((CX-ACT) )##2))/02.8(8(1)-A(1))) 
V=V4tVV 

CONTINUE 

IF(NOF)10,10,75 
75 CO 9 N=1,NOF 

IF (X+EPSX-XFO(N)) 9,958 
VV=FY(N) 
V=V4VV | 

9 CONTINUE 

10 IF(NOR)12,13,105 = 
105 DO 12 M=1,NOR 

LF (X+EPSX-XR(M))12,12,11 
11 VV= RY(M) 

V=V4+VV 
12 CONTINUE 
13 CONTINUE 

RETURN 

END 
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SUBROUTINE MOMENT (X, AM) 

DIMENSION A(25)28(25) ,QA(25),QB(25) »XFG(100) »FY(100) ,XR(1C)2RY(10) 
¥2XM(5C) ¢AMZ(50) »2XM(10),RAMZ(10) ,AE(25),8E(25),EA(25),EB( 25) 
29GA(25),G68( 25) ,AG (25) ,8G(25) 

COMMON A, BsQAsQB, XFOsFY2XRyRV ep XMy AMZ, RXMeRAMZ,AEZBE, EA, EB,NCL» 

1 NOR, NOFsNOMsNORM y NOI» NOK,» GA, GB, AGe BG 
EPSX=10.##(-8 ) 
AM=0. 
IF(NOL)8,8,2 
2 C00 7 fT=1,NOL 

IF(X-A(1))7,3,3 
3 IF(X+EPSX-B(1I))5, 554 
& AMM=(4+QA(1)/2.)#( (XACT) ) ##2)4+(° (CBU T)-QACT)) # CUX-ACT) ) # #3) )/(6 2% ( 
TBC I)—-ACT)))—(OB(I)/2.)8#0(X-BU I) )##2)— (QBCT)-CACT) )#((X-B (U1) )##3) 
2/(6.* (BUI )-ACTI))) 

GO TO 6 
S AMM=(4+QA(1)/2.)#( (XACT) )# #2) 40 (QB01) -QACT) )# ((X-A (CT) )##3))/(06.% 
TCBCTI-ACTI))) 
AM=AM+AMM 
CONTINUE 
IF(NOF) 12,1299 
CO 11 N=1,NOF 
- TF (X+EPSX-XFO(N)) 11,117,710 
1O AMM= FY(N)#*(X~XFO(N)) 

AM=AM+AMM 
11 CONTINUE 
12 TFCNOR) 16916213 
13 00 15 M=1,NOR 

IF(X+EPSX-XR(M))159175, 14 
Th AMM= RY(M)#(X-XR(M)) 

AM=AM+AMM 
15 CONTINUE 
16 IF(NOM)28,28,17 
17 CO 19 K=1,NOM 

IF (X+EPSX-XM(K) 119919, 18 
18 AMM=ANZ(K) 

AM=AM+AMM 
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CONTINUE 

IF (NORM) 32,32,29 

CO 31 KA=1,NORM 

IF CX+EPSX-RXM(KA) )31431,30 
AMM=RAMZ(KA) 

AM=AM+AMM 

CONTINUE 

CONTINUE 

RETURA 

ENO 
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SUBROUTINE SLOPE( X,Y1) 
DIMENSION A(25),8 (25),Q0A(25) ,QB(25) »XFC( 100) »FY(100) ,XR(10),RY(10) 
15XM(5C) AMZ (50) ,RXM( 10), RAMZ(10) »AE(25) »BE(25), EA 25) EB (25) 
21GA(25),GB(25),AG (25) »BG(25) 
COMMON A,B, QA, QBy XFO,FYsXRyRV_y XMyAMZ,RXMyRAMZ, AE SBE, EA, EB, NOL 
1 NOR,NOF,NOM,NORM sNOI »NOK,GA,GB,AG, BG 
Y1=0. 
BO 37 {A=1,NOI 
IF (X-AE(1IA))17317,2 
2 IF(NOL)5,5,3 
300 4 I=1,NOL 
20 EYY = EA(IA)*((QA(1)/2.) #ENTON(X,AE(IA) ACT) 92) #0 (CBC I)-QAC1) )/ 
1(6.8( BC 1)-ACT)))) #(ENTON(X,AE( TA) ACT), 3)-ENTON(X, AE( IA) »B(1) 23) ) 
2 -(QB(1)/2.)#ENTON( X,AE(TA),B(1),2)) +((EBCTA)-EA(TA)) /(BE( IA)~ 
BZ AE(IA)))#((OA(CT)/2.) #ENTIN(X,AEC IA) A(T) 92) #0(GBII)-CACI)) / 
U(6.#( BC T)—-ACI)))) *(ENTINGX,AEC IA) A(T)» 3)-ENTIN(X, AE( IA) B(1) 92) ) 
S -(QB(1)/2.)#ENTIN(X, AE( IA) B( 1) 92)) 
YI=Y1I+EYY | 
21 EYY= -EB(IA)#((QA(1)/2.) #ENTON(X,BE(IA) ,AC1),2) #0 (QB(I)-CACI) )/ 
1(60#(B(I)-A(T)))) #(ENTON( X, BEC IA) yA( 1) 43 )-ENTON(X, BE( IA), B(1),2)) 
2 -(QB(1)/2.) #ENTON(X,BE( IA) »6(1),2)) -(CEBCIA)-EA(IA))/(BE(IA)- 
3 AE(IA))) #( (QA(T)/2.) #ENTIN(X,BEC IA), ACI) ,2) +((CBII)-QA(I)) / 
U(6e#( BCT AC I) )))# (ENTIN( Xs BEC TA) pACT) 3) -ENTINGX, BE( IA) yB( 1)» 3)) 
S -(QB(1)/2.)*#ENTIN(X, BEC IA) »B(1),2)) 
YI=Y1I+EYY 
4u CONTINUE 
S IF(NOF)&@,8,6 
‘6 OO 7 N=1,NOF 
22 EYY= EA(IA) #FY(N) #ENTON(X,AE( IA), XFO(N) 21) -EBC IA) ®&FY(N) SENTON(X, 
1BE(IA),XFO(N),1) + (CEBCIA)-EACIA)) /(BECIA)-AE(IA) )) #FY(N)# 
2(ENTIN(X,AE( IA) »XFO(N) 21) — ENTIN(X,BE( IA), XFO(N),1)) 
YI=YI+EYY 
7 CONTIAUE 
8 IF(NOR)11,11;9 
9 00 10 M=1,NOR 
23 EYY= EA(IA)#RY(M) #ENTON(X,AE(TA),XR(M)21)-EBCIA) #RY(M)#® ENTON(X, 
 -TBEC IA) sXR(M) 91) # ( (EBC IA) -EACIA)) /(BEC(IA)-AE(IA)))#RY(M) « 
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2(ENTIN(X,AE( TA) XRIM) 9 1)— ENTIN(X, BEC IA) »XR(M)91)) 
YI=YI+EYY 
10 CONTINUE 
11 TFCNOM) 14, 14,12 
12 00 13. K=1,NOM | | 
26 EYY= EA(IA) #AMZ(K )#ENTON (Xs AE CTA) »XM(K),0) ~EB( IA) *AMZ(K) #ENTON(X, 
TBE( TA) »XM(K),0) + (CEB(ITA)-EA(TA))/(BE(TA)7AE(TA)) )® AMZ(K) 
2 (CENTINGXsAEC TA), XM(K),0)— ENTINEOXSBE( TA) »XM(K) 20) ) 
YI=VI+tEVY | 
13 CONTINUE 
14 ITFCNORM) 17,17%15 
15 DO 16 KA= 1,NORM | 
25 EYY=EA( 1A) #RAMZ(KA) *ENTON(X,AE( TA) »RXM(KA) »0)-EB (IA) *#RAMZ(KA)# 
1 ENTON(X,BECTA) SRXM(KA)20O) + (CEBC IA) -EACIADY/(BECIAD-AE(IA))) & 
2RAMZ(KA)#® CENT IN(X » AE (IA) »RXM(KA),0)-" ENTIN(X BEC IA) »RXM(KA) 70) ) 
YI=YI+EYY : 
16 CONTINUE ~ 
17 CONTIAUE 
IF (NGK)595, 595,500 
500 00 59G_ 18=1,NOK | 
IF (X -AG(1E))590,596,509 
509 IF (NCL)530,530,510 
510 06 52C 1=1,NOL 
26 EYY = GA(IB)#(QAC1)#FTON(X,AG(IB) pACI) 51) + (CQBUL) -QA(I) ) / 
1 (2.8 (BC L)-ACT))) ) #CETONOX,AG (IB) pAC I), 2)-FTON(X »AGCIB)  B( 1) 92)) 
2 -QB(1)#FTON(X,AG( IB) »B(1).1)). # ((GB¢IB)-GACIB) )/ (BG IB)-AG(IB))) 
3 #(QACL)®FTIN(X,AG(IB) ACT) 91) + CCQBCII-QACI)I/(20#(B(T)-AC1)))) 
bs (FTIN(X,AGCIB), ACI) 92) -FTIN(X,AG(IB) eB(1),2)) - QBCI)* FTIN 
S (X,AG(1B),B(1) 92)) 
Yl =¥1 -EYY | 
27 EYY =-GB(IB)#(QACI)#FTON(X,BGCIB),ACT),1) + ((QB(L) -CACI)) / 
1 (2.8 (BCL)-ACT))) )#€ FTON(X,BGCIB) sAC 1), 2)-FTON(X, BGC IB) ¢B(1),2)) 
.2-QB(1)* FTON(X,BG( 1B) ,B( 1) 91) )— ((GBCIB)-GAC(IB)) /(BG(IB)-AG(IB))) 
Ze(QACL)# FTIN(XsBGCIB) ACI) 91) (CQBCII-QACI)I/(2.#(B(1)-A(I)))) 
ye ( FTIN(X,BG(IB),A(I)92)— FTINGX,BG(1B),8(1)22)) - QBCI)# FTIN 
5 (X,BG(18),B(1)52)) | 
YI =Y¥1 -EYY 
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520 CONTINUE 
530 IF(NOF)56C,5602S40 
S40 DO 55C N=1,NOF 
28 EYY=GA(IB)#FY(N)* FTON(X,AG(IB),XFO(N),0) -GB(IB)#FY(N)® FTON(X, 
1BG(18),XFO(N) ,0) + ((GB(IB)-GA(IB)) / (BG(IB)-AG(IB))) * FY(N)# 
2( FTIN(X,AG(IB) »XFO(N) 20) - FTIN(X,8G( 18) »XFO(N),C)) 
Y1 =Y1! -EYY 
550 CONTINUE 
560 IF(NOR)590,590,570 
S70 DO 58C #=13,NOR 
29 EYY =GA(IB)#RY(M)# FTON(X,AG(IB),XR(M),0) -GB(I@)#RY(M)® FION(X, 
1BG(1B)»XR(M),0) + ((GB(IB)-GA(IB)) / (BG(IB8)-AG(1B))) * RY(M) # 
2( FTIN(X,AG(IB) »XR(M),0) -— FTIN(X,8G( IB) »XR(M),0C) ) 
Yl =Y1 -EYY 
5@0 CONTINUE 
590 CONTINUE 
595 CONTINUE 
RETURN 
END 
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' SUBROUTINE DEFLECT(X,Y) 

DIMENSION A(25) 5B (25) 4QA(25) yCB( 25) 5XFO( 100) »FY(100),XR(10)9RY(10) 
19XM(50)9AMZ(50) 9&2 XM(10),RAMZ( 10) yAE(25) 9 BE (25) yEA( 25), EB( 25) 
2+GA(25),GB(25)4AG (25) yBG (25) 

COMMON Ay By QA,QBy XFOsFY¢XRyRV yXMy AMZ» RXMy RAMZ » AEyBEyEAyEByNOLy 
1 NOR,NOF,NOM, NORM+ NOI »NOKy GA, GB, AGe BG 

T=0. | 

DO 17 IA=1,NOI 

IF (X—AE(IA})17,) 1742 

2 IF(NOL)5S,5,3 
300 4& I=1,NOL 
20 EYY = EAC IA)#((QA(1)/2.)*#DI TONG Xs AEC TA) ACT) 42) #¢(QB(I)-QA(1))/ 
1(60#(B(1)-A(T)) 2) @(DITON(X AEC TAD 9ACT) 9 3)-DITON(X, AE(IA) »B(1)93)) 
2 -(QBC1)/2.)#DITON(XsAE( TAI, BU 1)92)) + (EBC IAI-EAC IA) )/(BE(IA)- 
3 AE(IA)))#((QAC I) /2.)*#DITIN(X AEC IA)» ACT) 92) #((QB(T)-QA(I)) / 
&(6e#(B(T)—-ACI)))) #(CDITIN(XsAE(IA) pA(L) 93 )-DITING Xp AEC IA) »B(1)_3)) 
5-(Q8(1)/2.)#DITIN(X,AE(IA),B(I) 20) 
Y=Y+E YY 
21 EYY= -EBCIA)#((QA(T)/20)*DITON(XsBE( TA) A(T) 52) +( (QB(T)-QACT))/ 
1(60#(BCIT)-A(T)))) *# (GI TON(X,BE( TA)» ACT) ¢3)-DITON(X, BE( TA) BT), 3)) 
2 -(QB(1)/2.)#DITON(X,BE( IA), B(1),2)) —(CEBCIA)-EA(TA))/( BE(IA)- 
BAEC IA)) )#((QA(T)/ 20) *DITIN(X,BEGIA) A(T) 92) #((QB(T)-QAC1)) 4 
4 (6e#(B(I)-ACT)))) @(DITIN(X, BEC IA) sA(T) »3)-DITIN( X, BEC IA) »B(1),3)) 
5 -(QB(1)/2.)*#DITIN(X,BE(IA),B(1),2)) 
Y=Y+EYY i. 
&u CONTINUE 
5 IF(NOF)8,8,6 
6 DO 7 N=1,NOF 
22 EYY= EA(IA) #FY(N) #DITON(X,AE( IA) »XFO(N), 1) -EB( IA) #FY(N) ®DITON(X,, 
1 BE(IA),XFO(N),1) + ((EBCIA)-EA(TA))/(BE(IA)-AE(IA))) #FY(N)# 
2(DITIN(X, AEC IA) sXFOC(N) 91) — DITIN(XsBE( IA), XFO(N)»1)) 
Y=Y+EYY | 
7 CONTINUE 
8 IF(NOR)11,1179 
9 DO 10 M=1,NOR 
23 EYY= EACTA)*#RY(M) #DITON( X,AE( TA) ,XR(M), 1)-EB( IA) #RY(M)® DITON(X, 
IBE(IA),XR(M),3) +((EB(IA)-EACIA))/( BEC TAI-AE(IA))) #RY(M) # 
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2(DITIN(X, AE(IA) sXR(M)91)— DITIN(X,BE( IA) yXRUM) 91)) 
Y=Y+EYY 
10 CONTINUE 
11 TF(NOM) 14,14, 12 
12 DO 13. K=1,NOM 
24 EYY=EA( IA) *AMZ(K) #DITON(X,AE(IA) 9 XM(K),0) =EB(IA)#AMZ(K)#DITON(Xs 
IBECIA),XM(K),0) + (CEBCIA)-EA(TA))/(BEC(IA)<AE(IA)) )#AMZ(K)# 
2 (DITIN(X,AE( IA), XM(K) 90)— DITIN(X,BE (IA) »XM(K) 20) ) 
Y=Y+EYY 
13 CONTINUE 
14% IF(NORM) 17,17, 15 
15 DO 16. KA= 1,NORM 
25 EYY=EA(IA)#RAMZ(KA) #DITON(X,AE(IA) »RXM(KA) ,0)-EB(IA) #RAMZ (KA) # 
1 DITON(X,BE(IA)sRXM(KA) 90) + ((EBCIA)-EAC(IA))/(BECIA)-AE(IA))) # 
2RAMZ(KA)# (DIT IN(X pAE(IA) »RXM(KA) O)— DITIN(X,BE(IA) »RXM(KA)C)) 
Y=Y+EYY : 
16 CONTINUE 
17 CONTINUE 
1F(NOK)595,595, 500 
500 DO S9C_ IB=1,NOK 
IF (X -AG(1B))590,590,509 
509 IF (NCL)530,53C15 10 
510 00 52C =1,NOL 
26 EYY=GA(IB)#(QA(I) # ENTON(X,AG(IB),A(T),1) # ((QB(I) -CA(T) ) / 
1 (2e"8(B(L)-A(1))) )#CENTON (X,AG(IB) pACT) 92 )-ENTON(X pAG( IB) 9 B( 1) 92)) 
2~QB(1)*#ENTON(X,AG (IB) 9B(1)21)) + ((GB(IB)-GA(IB))/ (BG(IB)-AG(IB))) 
Ze (QACI)*ENTIN(XAGCIB) ACT) 91) + (CQBCI)-CACI)I/(2-#(B(1)-ACT)))) 
be (ENTIN(X,AG(IB), ACT), 2)-ENTIN(X,AG(IB) pB(1),2)) - QB(I) *#ENTIN 
5 (X,AG(IB),B(1),2)) 
Y=Y-EYY 
27 EYY=-GB(IB)#(QA(I )#ENTON(X,BG(IB),ACI),1) #((QR(I)-QA(I))/ 
1(2e#(BUI)-A(1)))) #(ENTON( X,BG(IB) A(T) »2)-ENTON(X, BG (IB) pB(1)92)) 
2-QB(1)*ENTON(X,8G (IB) ,B( I) 21) )— ((GBCIB)-GA(IB))/( 8G(IB)-AC(IB))) 
Ze (QACI)#ENTIN(X,BG(IB)sACT), 1+ ((QBCI)-QA(T))/(2. #(B(T)“A(T)))) 
he (ENTIN(X,pBG( IB), A(T) »2)-ENTIN(X? BG(1B),8(1),2)) - QB(1T)# ENTIN 
5 (X,BG(1B) ,B( I) 92)) 
Y=Y-EYY | 
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520 CONTINUE 
530 IF(NOF)560,5605540 
540 CO 550 N=1;NOF 
28 EYY=GA(IB) *#FY(N)* ENTON(X,AG(IB) »XFO(N) 20) -GB(IB)#FY(N) #ENTON(X, 
1BG(18),XFO(N) »0) + ((GB(IB)<GA(IB)) / (BG(IB)-AG(I8))) * FY(N)« 
2 (ENTIN(X,AGCIB)sXFO(N),0) — ENTIN(X,BGC(IB) »XFOC(N),C)) 
Y=Y-EYY 
550 CONTINUE 
560 IF(NOR) 590,590,570 
570 DO 580 M=1,NOR 
29 EYY =GA(IB)#RY(M) *ENTON(X,AG(IB) »XR(M),0) -GB(IB)#RY(M) *ENTON( Xy 
1BG(IB),XR(M),0) + ((GBCIB)-GA(IB)) / (BG(IB)-AG(IB))) * RY(M) * 
2(ENTIN(XsAG(IB)sXROM),0) — ENTIN(X,8G( 1B) yXRO(M) 0) ) 
Y=Y-EYY 
580 CONTINUE 
590 CONTINUE 
595 CONTINUE 
RETURN 
END 
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FUNCTION FTON(X yA yBgN) 
FTON = 0. 

HeeX— A013: 351 

TF UX-01303% 2 
FTON=(X-B) ¥«N 

RETURN 

END 


FUNCTION FTIN (X,A,B,N) 


FTIN =C. 


IF (X-A)3,3,1 

IF (X-B) 3932 

FTIN = (X-B)##(N#17) + (B-A)#((X-B) ##N) 
RETURN 

END 


\ 
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FUNCTION ENTON (X,A,B,N) 


FN=N 

ENTON = O. 

IF(X-A)5,5, 1 
1 IFCX-8)5,5,2 : 
2 IF(A-8)3,3,4 
3 ENTON = (X=-B)##(N41)/(FN41.) 

GO TC 5 

ENTON = (X-B)##(NtT)/(FN4+12) = (A=-8) ##(N41)/(0FN41.) 
5 RETURN 

END 


FUNCTICGCN ENTIN( X,A,8,N) 

FN=N 

ENTIN = O. 
IF (X-4)5,5,1 
IF (X-B)5,5,2 

IF (A-8)3,3,4 

3 ENTIN = (X-B)##( N42) /(FNt2.) + (BrA)#(X-B)##(N41)/(0F N41.) 
. GO TC 5 

GY ENTIN = (X=-B)##(N4#2)/(FN4t2.) + (BrA) #(X—B8) ##(N+1) SC FN4+1.) 

1 =-(A-B) ##(N#2)/(F N42.) —( (BA) #(A-B) ##(N4+1))/ (FN41.) 
5 RETURN Z 
ENO 


ie 
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FUNCTION DITON(X,A,B,N) 
FN=N 
OITON =0. 
IF(X-A)5,5,1 
IF (X-2)5,53;2 
2 IF(A-8)3,3,4 


3 OITON = (X-B)##(N+2)/(°00FN41.)#(FN+2.)) 
GO TO 5 
4 DITON = (X-B) ##(N42)/0(FN+1.)4(0FN4t2.0)) —(X#(A-8) ##(N41))/0FN+1.) 


1-(A-B) #¥(Nt2)/((F NT.) #(FN42.2)) + (A®(A-B) ## (N41) )/(CFNt+1.) 
5 RETURN 
END 


FUNCTION DITIN(X»Ae8,N) 
FN=N 
DITIN=0. 
IF(X-A) 5,5, 1 
1 IF(X-8)5,5,2 
2 IF(A-B)3,3,4 
3 DITIN =6X—B) ##(Nt3)/( (FNt2.) 40 FNt3.)) + ((B-A)#(X—B) HA(N42))/ 
1((CFN+1.)#(FN4+2.)) 
GO TO 5 7 
W DITIN = (X-B) ##(N43)/0(FN42.) #(FN4+3.)) + ((B-A)#(X-B) eH (N42) )/ 
10 CEN+ 1.) #(FN42.)) —(X#(A-B) ##(N+2)) JC EN42.)—-(X# (BA) #(A-B) ##(K41)) / 
2 (FN41.) — (AB) # #(N43)/0(0F N42.) #(FN432)) -— (BA) #(A-B)##(NH2) / 
3 CCFEN4+1.)@(FN42.)) + (A®(A-B) ##(N42))/(0FN42.) + (A#(B-A) #(A-B) 2% 
GW (N+1))/(0FN+1.) 
S RETURN 
END 
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SUBROUTINE GAUSS2(N,EP,A,X,KER) 


CIMENSION A22,23) »X(22) 


NPM=N+1 

DO 34 L=1,N 

KP #Q 

Z=0.C 

DO 12 K=L,N 

IFC Z~ABSF(A(K,L) DV 115125192 
Z=ABSF(AC(K,L)) 

KP=K 

CONTINUE 

IF(L-KP) 13,20, 20 


“00 14 J=L,NPM 


Z=A(LyJ) 

A(LyJ)=A(KP, J) 

A(KP,J)=Z 

IF (ABSF(A(L,L) )-EP)50,50, 30 
IF(L=N)31,40,40 ‘ 

LP 1=L+] 

DO 34 K=LPI1,N 
IF(A(KyL)) 32, 24, 32 
RATIO=A(KyL)/A(LyL) 

00 32 J=LP1,NPM 
A(K,J)=A(KyJ)—-RATIO#A(Ls J) 


CONTINUE oh 


DO 4&2 I=1,N 
IT=N+1~I1 

DO 42 J=l1,] 
JPN=J+4+N 

$=0.0 

IFC TIA“N)81,43,4%3 
TIPI=II+1 

DO 4&2 K=IIPI,N 
S=StAC II ,K)#X(K) 
X(II) =CACII,JPN) -S)I/ACTII,IT) 
KER=1 


RETURN 
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50 KER=2 
RETURN. 
END 
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